


clear
capture log close
set more off
log using growth_inequality.smcl, replace
global path ""
cd ""
global dta1 ""



*** recs 1993 data as an alternative robustness check
******* cpi 2
. clear

. use "bls_prices_2.dta"

* normalize cpi to equal one in base years

foreach num in 1993  {
	gen norm`num' = cpi if year==`num'
	sort year
	replace norm`num' = norm`num'[_n-1] if norm`num'==.
	gsort -year
	replace norm`num' = norm`num'[_n-1] if norm`num'==.
	gen cpi_`num' = cpi/norm`num'
	
}
* save
keep year cpi_*
sort year
save cpi_allyears, replace 

clear
use cps_00003
sort year serial pernum
merge m:1 year using cpi_allyears
*keep if _merge==3
drop _merge

* calculate real income

gen hhincome_1993 = hhincome/cpi_1993

** keep household head only
keep if relate==101
	
* Census division
gen division = 1 if region==11
replace division = 2 if region==12
replace division = 3 if region==21
replace division = 4 if region==22
replace division = 5 if region==31
replace division = 6 if region==32
replace division = 7 if region==33
replace division = 8 if region==41
replace division = 9 if region==42

save income_b1993.dta,replace


***0918
clear
use income_b1993.dta

forvalues i=1993/2020 {
use income_b1993.dta

* create real income groups

	gen hhfaminc_1993 = 1 if hhincome_1993<3000
	replace hhfaminc_1993 = 2 if hhincome_1993>=3000 & hhincome_1993<4000
	replace hhfaminc_1993 = 3 if hhincome_1993>=4000 & hhincome_1993<5000
	replace hhfaminc_1993 = 4 if hhincome_1993>=5000 & hhincome_1993<6000
	replace hhfaminc_1993 = 5 if hhincome_1993>=6000 & hhincome_1993<7500
	replace hhfaminc_1993 = 6 if hhincome_1993>=7500 & hhincome_1993<9000
	replace hhfaminc_1993 = 7 if hhincome_1993>=9000 & hhincome_1993<10000
	replace hhfaminc_1993 = 8 if hhincome_1993>=10000 & hhincome_1993<11000
	replace hhfaminc_1993 = 9 if hhincome_1993>=11000 & hhincome_1993<12500
	replace hhfaminc_1993 = 10 if hhincome_1993>=12500 & hhincome_1993<14000
	replace hhfaminc_1993 = 11 if hhincome_1993>=14000 & hhincome_1993<15000
	replace hhfaminc_1993 = 12 if hhincome_1993>=15000 & hhincome_1993<17500
	replace hhfaminc_1993 = 13 if hhincome_1993>=17500 & hhincome_1993<20000
	replace hhfaminc_1993 = 14 if hhincome_1993>=20000 & hhincome_1993<22500
	replace hhfaminc_1993 = 15 if hhincome_1993>=22500 & hhincome_1993<25000
	replace hhfaminc_1993 = 16 if hhincome_1993>=25000 & hhincome_1993<27500
	replace hhfaminc_1993 = 17 if hhincome_1993>=27500 & hhincome_1993<30000
	replace hhfaminc_1993 = 18 if hhincome_1993>=30000 & hhincome_1993<32500
	replace hhfaminc_1993 = 19  if hhincome_1993>=32500 & hhincome_1993<35000
	replace hhfaminc_1993 = 20 if hhincome_1993>=35000 & hhincome_1993<40000
	replace hhfaminc_1993 = 21 if hhincome_1993>=40000 & hhincome_1993<45000
	replace hhfaminc_1993 = 22 if hhincome_1993>=45000 & hhincome_1993<50000
	replace hhfaminc_1993 = 23 if hhincome_1993>=50000 & hhincome_1993<75000
	replace hhfaminc_1993 = 24 if hhincome_1993>=75000 & hhincome_1993<100000
	replace hhfaminc_1993 = 25 if hhincome_1993>=100000 & hhincome_1993!=.

keep if year==`i'

collapse (sum) asecfwt, by(hhfaminc_1993 division)
rename asecfwt hh_count

egen sumhh=total(hh_count)
gen weight=hh_count/sumhh
save `i'_93_division,replace
}



*** predicted kwh

clear
forval i=1993/2020  {
clear
use recs_1993

bysort division moneypy: egen sumhh=total(nweight)
bysort division moneypy: gen weight=nweight/sumhh
bysort division moneypy: egen meankwh=total(weight*kwh)

 collapse meankwh,by(moneypy division)

quietly drop if division==. 
rename moneypy hhfaminc_1993
sort hhfaminc_1993 division
quietly destring hhfaminc_1993,replace

quietly merge 1:1 hhfaminc_1993 division using  `i'_93_division

quietly drop if meankwh==.
quietly drop if hh_count==.

*sum kwh // this sum is not estimated/ estimated is by weights, this is just summary
quietly sum meankwh[aweight=weight]
quietly return list
dis `i',r(mean)

}

*****====standard errors====*****

clear
use recs_1993

sort moneypy division

egen hhgroup=group(moneypy division)

reg kwh i.hhgroup [weight=nweight], nocons //bootstrap seems not allow weights
estimates store results1
esttab results1, b(3) se(3) scalars(r2) star( * 0.10 ** 0.05 *** 0.01) 


matrix B=e(V)

local dim (`= rowsof(A)',`=colsof(A)') 

di "`dim'" 

forval i=1993/2020 {
quietly use `i'_93_division,clear
quietly mkmat weight,matrix(A)
quietly matrix SD=A'*B*A
quietly matrix list SD
quietly scalar z=SD[1,1]
di sqrt(z)
}

